\documentclass[twocolumn]{article}
%\documentclass[]{article}
\usepackage{amsmath,amsfonts,afterpage}
%\usepackage{showlabels}
\usepackage{lscape}
\usepackage[pdftex]{graphicx,color}
\newcommand{\normal}[2]{{\cal N}(#1,#2)} \newcommand{\La}{{\cal L}}
\newcommand{\nomf}{\tilde f} \newcommand{\COST}{\cal C}
\newcommand{\LL}{{\cal L}} \newcommand{\Prob}{\text{Prob}}
\newcommand{\field}[1]{\mathbb{#1}}
\newcommand\REAL{\field{R}}
\newcommand\Z{\field{Z}}
\newcommand\Polytope[1]{\field{P}_{#1}}
\newcommand\PolytopeN{\Polytope{N}}
\newcommand\PolytopeInf{\Polytope{\infty}}
\newcommand{\EV}[2]{\field{E}_{#1}\left[#2\right]}
\newcommand{\partialfixed}[3]{\left. \frac{\partial #1}{\partial
      #2}\right|_#3}
\newcommand{\partialtwo}[2]{ \frac{\partial #1}{\partial
      #2}}
\newcommand{\dhot}{=_{\text{dhot}}}
\newcommand{\argmin}{\operatorname*{argmin}}
\newcommand{\argmax}{\operatorname*{argmax}}
\newcommand{\set}[1]{\left\{ #1 \right\}}
\newcommand\inner[2]{\left<#1,#2\right>}
\newcommand\norm[1]{\left|#1\right|}
\newcommand\bv{\mathbf{v}}
\newcommand\bw{\mathbf{w}}
\newcommand\dum{\xi}
\newcommand\Ddum{d\dum}
\newcommand\lambdaPF{\lambda_{\rm{PF}}} % Peron Frobenius eigenvalue
\newcommand\rhoPF{\rho_{\rm{PF}}} % Peron Frobenius eigenvector

\newtheorem{conj}{Conjecture}

\input{eric.latex} % newcommands from sympy calculations

\title{Quadratic Approximation to First Order Integral Operator}

\author{Andrew M.\ Fraser}
\begin{document}
\maketitle
\begin{abstract}
  Based on \emph{Propagating Uncertainty About Gas EOS to Performance
    Bounds for an Ideal Gun} LA-UR-12-22731, I describe simulations
  and calculations on an approximation to an
  integral operator.\\
  \textbf{To do:}
  \begin{itemize}
  \item Revise simulation code to assume $\Delta_y = 1$, and
    numerically check Conjecture~\ref{conj}
  \item Make alternative plans if conjecture not true
  \item Refactor and submit for publication.
  \end{itemize}
\end{abstract}
\section{Introduction}
  \label{sec:introduction}

Here I consider a first order Markov process on 2-d states consisting
of function values $g(y)$ and derivatives $h(y) \equiv \left. \frac{d
    g }{d t} \right|_{t=y}$.  I consider sequential samples at $y_0$
and $y_1$, and I suppose that the upper bound is $u$ and the lower
bound is $l=-u$.

In $(f,x)$ coordinates, I describe the nominal function, $\tilde f$ by
\begin{align*}
  x_0 &= 1\\
  f_0 &= f(x_0) \\
  \tilde f(x) &= \frac{\tilde f_0}{x^3}.
\end{align*}
The coordinates $g$ and $y$ are functions of coordinates $f$ and $x$ with
\begin{align*}
  x(y) &= e^y \\
  f(g) &= \frac{\tilde f_0 e^{g}}{x^3} \\
  g(f) &= \log(f) + 3y - \log(\tilde f_0).
\end{align*}
The straight line in $(g,y)$, $g(y) = 0$, corresponds to
\begin{equation*}
  \tilde f(x) = \frac{\tilde f_0}{x^3}.
\end{equation*}

In the $(f,x)$ coordinates a line that is tangent to $\tilde f(x)$ at
$x=1$ satisfies
\begin{equation*}
  f(x) = f_0(4-3x),
\end{equation*}
which in the $(g,y)$ coordinates is
\begin{equation*}
  g(y) = \log \left( 4 - 3e^y \right) + 3y.
\end{equation*}
At $y=0$, that $g$ has the following derivatives
\begin{align*}
  g(0) & = g_0 \\
  g'(y) &= 3 - \frac{3e^y}{4-3e^y} \\
  g'(0) &= 0 \\
  g'' &= -\frac{\frac{4}{3}e^y}{\left( e^y - \frac{4}{3} \right)^2}\\
  g''(0) &= -12 \\
  g''' &= -\frac{4\left( \frac{4}{3} + e^y \right) e^y}{
    3\left( \frac{4}{3} - e^y \right) ^3} \\
  g'''(0) &= -84.
\end{align*}
Thus the second order Taylor series approximation to the function in
$(g,y)$ that corresponds to a tangent to $\frac{\tilde f_0
  e^{g_0}}{x^3}$ at $x_0 = 1$ is
\begin{equation}
  \label{eq:tangent}
  g(y) = g_0 - 6(y-b)^2.
\end{equation}
I will use the Taylor series approximation \eqref{eq:tangent} to
tangents in $(f,x)$ from now on.  Figure~\ref{fig:taylor} illustrates
that approximation.

\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{taylor.pdf}}  
    \caption{Second order Taylor series approximation to $g(y)$
      compared to actual $g(y)$ when $f = e^g = \alpha - \beta x$,
      $y=\log(x)$, and $f$ is tangent to $\frac{1}{x^3}$ at $x=1$.
      The indicated bounds, $\pm 0.0005 = 5\times10^{-4}$ correspond
      to a factor of $1.0005$ which is one tenth the value quoted by
      Shaw.}
  \label{fig:taylor}
\end{figure}
%\end{document}
Now I describe quadratic approximations to constraints on the values
of $g_1 \equiv g(y_1)$ and $h_1 \equiv h(y_1)$ that are determined by
$g_0 \equiv g(y_0)$ and $h_0 \equiv h(y_0)$.  I use $U_g$, $L_g$,
$U_h$, and $L_h$ to denote (upper/lower) bounds on ($g_1/h_1$)
respectively.  A trajectory $g$ that goes through $g_0$ at $y_0$ must
lie between the two lines $U_g(g_0,y_0,y)$ and $\bar U_g(g_0,y_0,y)$
that in $(f,x)$ coordinates go through $f_0$ and are tangent to the
upper bound $u$ to the right and left of $y_0$ respectively.  If $h_0$
is the derivative of $g$ at $y_0$, the value of $g$ at $y_1$, ie,
$g_1$ must lie above $L_g(g_0,h_0,y_0, y_1)$, the image in $(g,y)$ of
the tangent to $f$ at $x_0$ in $(f,x)$.  Thus
\begin{equation}
  \label{eq:boundsA}
  U_g(g_0,y_0, y_1) \geq g_1 \geq L_g(g_0,h_0,y_0, y_1).
\end{equation}
The following calculations fit quadratic approximations to both $U_g$
and $L_g$.
\newcommand{\Rad}[1]{\sqrt{24\left(u#1\right)}}
\newcommand{\UgQ}{ g_0 + \Delta_y\Rad{-g_0} - 6 \Delta_y^2}
\begin{align}
  U_g(y) &= u - 6(y-b)^2 &&\text{Premise} \nonumber \\
  U_g(y_0) &= g_0 &&\text{Premise} \nonumber \\
  \label{eq:premise}
  U'_g(y_0) &\geq 0 &&\text{Premise} \\
  g_0 &= u - 6(y_0-b)^2 \nonumber \\
  \label{eq:root}
  b &= y_0 \pm \sqrt{\frac{u-g_0}{6}} \\
  \label{eq:dugb}
  \frac{d U_g}{d y} &= -12(y-b) && \text{Positive root in \eqref{eq:root}} \\
  \label{eq:dug}
  \left. \frac{d U_g}{d y} \right|_{y_0} &= 12\sqrt{\frac{u-g_0}{6}}
\end{align}
\begin{equation*}
  U_g(y) = g_0 + (y-y_0)\Rad{-g_0} - 6 (y-y_0)^2
\end{equation*}

Thus
\begin{equation}
  \label{eq:ug}
  U_g(y_1) = \min \left( u,\UgQ \right)
\end{equation}
Similarly, for the lower bound:
\newcommand{\LgQ}{g_0 + h_0\Delta_y - 6\Delta_y^2}
\begin{align}
  L_g(y) &= a - 6(y-b)^2 &&\text{Premise} \nonumber \\
  L'_g(y) &= -12(y-b) \nonumber \\
  g_0 &= a - 6(y_0-b)^2 &&\text{Premise} \nonumber \\
  h_0 &= -12(y_0-b) &&\text{Premise} \nonumber \\
  b &= y_0 + \frac{h_0}{12} \nonumber \\
  a &= g_0 + \frac{h^2_0}{24} \\
  \label{eq:lg}
  L_g(y_1) &= \max \left( -u,~\LgQ \right)
\end{align}
Figure~\ref{fig:boundsA} illustrates the bounds
\eqref{eq:boundsA}.  In Fig.~\ref{fig:boundsB}, I have simply reduced
the extent in $y$.

\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{bounds_04.pdf}}  
  \caption{Bounds on $g_1$ given $g_0$ and $h_0$}
  \label{fig:boundsA}
\end{figure}
\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{bounds_005.pdf}}  
  \caption{Bounds on $g_1$ given $g_0$ and $h_0$}
  \label{fig:boundsB}
\end{figure}

Given $g_0$ and a $g_1$ that satisfies the constraints
\eqref{eq:boundsA}, the constraints on $h_1$ are that it must be less
than the slope $\left. \frac{d U_g(g_1,y_1,y)}{dy} \right|_{y=y_1}$
and it must be greater than the slope at $y_1$ of the image in $(g,y)$
of the line in $(f,x)$ that connects the images of $(g_0,y_0)$ and
$(g_1,y_1)$.  The following calculation finds the quadratic
approximation for that curve:
\begin{align*}
  g(y) &= a -6(b-y)^2 \\
  g_0 &= a -6(b-y_0)^2 \\
  g_1 &= a -6(b-y_1)^2 \\
  a &= g_0 + 6(b-y_0)^2 \\
  &=  g_1 + 6(b-y_1)^2 \\
  g_0 -g_1 &= 12b(y_0-y_1) + 6(y_1-y_0)(y_1+y_0) \\
  b &= \frac{g_1-g_0 + 6 (y_1^2-y_0^2)}{12(y_1-y_0)} \\
  &= \frac{g_1 -
    g_0}{12 (y_1 - y_0)} + \frac{y_1 + y_0}{2}\\
  a &= \frac{\Delta_g^2}{24 \Delta y ^2} + \frac{g_1 + g_2}{2} +
  \frac{ 3 \Delta_y^2}{2}.
\end{align*}
At $y_1$, the derivative of $a-6(y-b)^2$ is
\newcommand{\LhQ}{\frac{\Delta_g}{\Delta_y} - 6 \Delta_y}
\begin{align}
  L_h(g_1, g_0, y_0, y_1) &= -12(y_1-b) \nonumber \\
  &= \frac{g_1 - g_0 -6(y_1 - y_0)^2} {y_1-y_0} \nonumber \\
  \label{eq:lower_h}
  & \equiv \LhQ
\end{align}
Following \eqref{eq:dug}, the upper bound is
\newcommand{\Uh}{\Rad{-g_1}}
\begin{equation}
  \label{eq:upper_h}
  U_h(g_1) = \Uh,
\end{equation}
and\footnote{The \emph{maximum} function and the term $-\Uh$ treat the
case of $U_g(y)$ being tangent to $u$ between $y_0$ and $y_1$.}
\begin{equation}
  \label{eq:UL_h}
  \Uh \geq h_1 \geq \text{max}\left( \LhQ,~ -\Uh \right).
\end{equation}

Figure~\ref{fig:boundsC} illustrates those constraints.

\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{bounds_dg.pdf}}  
    \caption{Curves $U_h(g_1)$ and $L_h(g_1, g_0)$ whose tangents give
      bounds on $h_1$ given $g_0$ and $g_1$}
  \label{fig:boundsC}
\end{figure}

\section{Scaling and Pie Slices}
\label{sec:scaling}

While both $\Delta_y$ and $u$ appear in \eqref{eq:ug}, \eqref{eq:lg}
and \eqref{eq:UL_h}, the ratio $d = \frac{u}{\Delta_y^2}$ is the
number that characterizes the problem.  The change of coordinates
\begin{subequations}
  \label{eq:rescale}
  \begin{align}
    \tilde g(g) &= \frac{g}{\Delta_y^2} \\
    \tilde h(h) &= \frac{h}{\Delta_y},
  \end{align}
\end{subequations}
is equivalent to substituting $d \rightarrow u$ and
$\Delta_y \rightarrow 1$ in the original problem.  In the new
coordinates, one finds that the structure of allowed sequential pairs
near $\tilde g=d$ depends only on the difference $d-\tilde g$ and not
on $d$ itself.  Without loss of generality, in this section, I assume
$\Delta_y=1$ and let $d$ define the problem.  I will use $h$ and $g$
as variables rather than $\tilde h$ and $\tilde g$.

The domain of the problem is the set of points
\begin{equation}
  \label{eq:domain}
  D = \left\{ z = (h,g) :
    \begin{array}{l}
      -d \, \leq g \, \leq d \\
      h^2\, \leq \, 24(d-g)
    \end{array}
    \right\},
\end{equation}
and the set of possible successors of each point $z_0 \in D$ is a
\emph{pie slice} $P(z_0)$ defined by the intersection of a wedge and
$D$.  In Figure~\ref{fig:eric}, I've labeled the points\footnote{I
  derived these expressions and the formulas for the ellipses in
  Figure~\ref{fig:eric} using the sympy symbolic calculation package.}
\begin{align*}
  z_1 &= \left( \hb, \gb \right) \\
  z_2 &= \left( \hc, \gc \right) \\
  z_3 &= \left( \hd, \gd \right)
\end{align*}
where $z_1$ is the wedge apex and $z_2$/$z_3$ is the intersection of
the upper/lower edge of the wedge with the boundary
respectively\footnote{Note that since the diagonal edge of each wedge
  goes through $(-6,g_0)$, $z_2$ depends on $g_0$ but not $h_0$.}.  In
summary, $z_1$ is an allowed successor of $z_0$ if both
\begin{subequations}
  \label{eq:symallowed}
  \begin{align}
    g_0 + h_0 -6 &\leq g_1 \leq g_0 + \sqrt{24(d-g_0)} - 6 \\
    g_1 - g_0 &\leq h_1 \leq \sqrt{24(d-g_1)}
  \end{align}
\end{subequations}

In Figure~\ref{fig:eric}, I use the ellipses that are the level sets
\begin{equation}
  \label{eq:moments}
  \left\{ z : (z-\mu)^T \Sigma^{-1} (z-\mu) = 2 \right\}
\end{equation}
for the image sets $P(z_0)$ of four initial points $z_0$ to illustrate
the mean $\mu$ and covariance $\Sigma$ of each pie slice.

\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{eric.pdf}}  
    \caption{Pie slice images of four points and ellipses that
      illustrate the means and covariances of those sets (See
      Eqn.~\eqref{eq:moments}).}
  \label{fig:eric}
\end{figure}

The following simple formula describes iterations of the map $z_1 =
F(z_0)$:
\begin{align*}
  F(h,g) &\equiv (h-12, g + h -6) \\
  F^n &\equiv (h_n, g_n) \\
  h_n &= h - 12 n \\
  g_{n+1} &= g_n + h - 12 n -6 \\
  g_{n+1} - g_n &= h - 12 n -6 \\
  g_{n+1} &= g + (n+1) (h-6) - 12\sum_{k=1}^n k \\
  &= g + (n+1)(h - 6(n+1)) \\
  g_n &= g + n(h - 6 n).
\end{align*}
In Figure~\ref{fig:f_n}, the boundary and the set
\begin{equation*}
S = \left\{ F^{\frac{n}{5}} (0,d-2.5) : n
    \in \left\{ -20, -19, \ldots , 20 \right\} \right\}  
\end{equation*}
appear.  Notice that each element of $S$ is below the boundary by the
same distance, namely 2.5.
\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{f_n.pdf}}  
    \caption{Illustration of the action of $F^n$ with noninteger
      values of $n$ and the action of $F$ on a set of points on the
      center, $h=0$.}
  \label{fig:f_n}
\end{figure}

The following calculation proves that $F$ preserves distance to the
boundary in the $g$ direction.  For a given $h$ the value of $g$ at
the boundary is
\begin{equation*}
  B(h) = d - \frac{h^2}{24}.
\end{equation*}
Now
\begin{equation*}
  F(h,B(h)-c) = \left( h-12, d - \frac{h^2}{24} - 6 -c \right)
\end{equation*}
and since
\begin{equation*}
  B(h-12) = d - \frac{h^2}{24} - 6,
\end{equation*}
The distance to the boundary in the $g$ direction, namely $c$, is
preserved.

\section{The Volume of Allowed Paths}
\label{sec:volume}

Here, I define some notation and quote an expression for joint
probability from Shannon.%
\newcommand{\tuple}[2]{z_{#1}^{#2}}%
\newcommand{\allowed}{\alpha}%
\newcommand{\preimage}{\tau}%
\newcommand{\scale}[1]{\psi_{#1}}%
\newcommand{\stationary}{\mu}%
\newcommand{\pair}{\beta}%
\newcommand{\ptuple}[1]{\gamma_{#1}}%
\newcommand{\Skip}[1]{\delta_{#1}}%
\newcommand{\PFop}{\Phi} %Peron Frobenius operator
\newcommand{\pfop}{\phi} %Continuum limit
\newcommand{\lefunc}{\rho_L}%
\newcommand{\refunc}{\rho_R}%
\newcommand{\eval}{\lambda}
\newcommand{\suparg}[1]{^{[#1]}}
\newcommand{\supn}{\suparg{n}}
\begin{description}
\item[$z\in D$] A point in the domain defined in \eqref{eq:domain}
\item[$\tuple{n}{m}$] A sequence $(z_n, z_{n+1}, \ldots, z_m)$
\item[$\allowed$] Indicator function
  \begin{equation*}
    \allowed(z_a,z_b) =
    \begin{cases}
      1 & \text{if }z_a \mapsto z_b
      \text{ is allowed (See \eqref{eq:symallowed})}\\
      0 & \text{otherwise}.
    \end{cases}
  \end{equation*}
\item[$\PFop$] The Peron Frobenius operator which acts on a
  function $f:D\mapsto \REAL$ by
  \begin{subequations}
    \label{eq:PFop}
    \begin{align}
      \left[f \cdot \PFop \right] (z) &\equiv \int f(z') \allowed(z',z)
      dz' \\
      \left[ \PFop \cdot f \right] (z) &\equiv \int \allowed(z,z')
      f(z') dz'.
    \end{align}
  \end{subequations}
  See \eqref{eq:domain} for the definition of $D$. 
\item[$\eval$, $\refunc$] The eigenvalue and the right eigenfunction
  of $\PFop$ defined by
  \begin{equation*}
    [\PFop \cdot \refunc] (z) = \eval \refunc(z)
  \end{equation*}
  I assume that the right and left eigenfunctions $\refunc$, $\lefunc$
  are normalized so that the joint probability function defined in
  \eqref{eq:pairdef} integrates to one.
\item[$\pair$] The joint probability density function for sequential
  pairs of points.  Theorem 8 of Shannon's 1948 paper says
  \begin{equation}
    \label{eq:pairdef}
    \pair(z_a,z_b) = \lefunc(z_a) \allowed(z_a,z_b) \refunc(z_b)
  \end{equation}
\item[$\stationary$] The stationary probability density
  \begin{align}
    \stationary(z_a) &= \int \pair(z_a,z) dz = \int \pair(z,z_a) dz
                       \nonumber \\
    \label{eq:stationary}
                     &= \eval \lefunc(z_a) \refunc(z_a).
  \end{align}
\item[$\ptuple{n}$] The joint probability density function for
  sequences of $n$ points
  \begin{equation}
    \label{eq:ptuple}
    \ptuple{n}(\tuple{1}{n}) = \frac{\lefunc(z_1)
      \refunc(z_n)}{\eval^{n-2}} \prod_{i=1}^{n-1} \allowed(z_i,z_{i+1})
  \end{equation}
\item[$\Skip{n}$] The joint density of pairs separated by $n$ skipped
  values, ie,
  \begin{align}
    \Skip{n}(z_a,z_b) &= \frac{\lefunc(z_a) \refunc(z_b)}{\eval^n}
                        \times \nonumber \\
    \label{eq:skip}
                      & \int \allowed(z_a,z_1) \prod_{i=1}^{n-1}
                        \allowed(z_i,z_{i+1}) \allowed(z_n,z_b) d
                        \tuple{1}{n}.
  \end{align}
  Note that the integral in \eqref{eq:skip} can be written as
  \begin{align}
    V_n(z_a,z_b) &\equiv \int \left[ \allowed(z_a,\cdot) \cdot
      \PFop^{n-1} \right] (z)~ \allowed(z,z_b)~ dz \nonumber \\
    \label{eq:volume}
    &= \left[ \allowed(z_a,\cdot) \cdot \PFop^{n} \right] (z_b),
  \end{align}
  and that it is the volume in $\Z^n \subset \REAL^{2n}$ of the set of
  allowed paths from $z_a$ to $z_b$ in $n$ steps.
\item[$\scale{n}$] A rescaling obtained by dividing $\Delta_y$ by $n$,
  ie, replacing $d$ with $n^2d$.
\item[$z\supn, \allowed\supn,\cdots $] Corresponding objects for
  rescaled problem, eg, if $z=(h,g)$ then
  \begin{equation*}
    z\supn = \scale{n}(z) = (nh, n^2g).
  \end{equation*}
  Note that $\allowed\supn$ is like $\allowed$ except that
  $\allowed\supn$ operates on a larger domain.
\end{description}

\subsection{Conjecture (I Don't Believe it Anymore)}
\label{sec:conjecture}

From \eqref{eq:skip} it seems that I need some kind of convergence of
$\PFop^n$.  I will define the required convergence in this section and
make a formal conjecture that it exists.

Let $\PFop\suparg{n,m}$ denote a map from functions on $D$ to
functions on $D$ defined so that starting with a function $u$, one
calculates $v = \PFop\suparg{n,m}u$ as follows
\begin{enumerate}
\item Map $u$ to higher resolution ($\Delta_y \rightarrow
\frac{\Delta_y}{n}$ or equivalently $d \rightarrow n^2 d$) with
$u' = \scale{n}(u)$
\item Apply the one step operator $\PFop$ appropriate for that scale
  $m$ times to yield $v' = \PFop^m(u')$
\item Map back to lower resolution with $v = \scale{\frac{1}{n}}(v')$
\end{enumerate}

Roughly, I want to define $\pfop(r)$ to be the continuum limit for
$y_2 - y_1 =r$ with
\begin{equation}
  \label{eq:phir}
  \pfop(r) = \lim_{n \rightarrow \infty} \PFop\suparg{n,rn}.
\end{equation}

Although I neglected to constrain the properties of $u$ in my
description of $\PFop\suparg{n,m}$, I only use the operator to
calculate volumes as in \eqref{eq:volume}.  So I simply define
\begin{equation}
  \label{eq:def}
  v_{n,m}(z_0,z) \equiv C \left[
    \allowed\supn \left( (z_0), \cdot \right)
    \cdot \PFop\suparg{n,m}
  \right](z),
\end{equation}
where $C$ is defined to ensure that
\begin{equation*}
  \int_D v_{n,m}(z_0,z) dz = 1.
\end{equation*}
For a particular $z_0$ the value of $v_{n,m}(z_0,\cdot)$ at $z$ is
proportional to the volume in $\REAL^m$ that corresponds to allowed
paths from $z_0$ to $z$.  As $m$ and $n$ get larger with a fixed
ratio, I hope that $v_{n,m}(z_0,\cdot)$ converges.  The following
conjecture is sufficient:
\begin{conj}
  \label{conj}
  \begin{align*}
    & \forall \epsilon > 0,~ \exists n_0 \text{ such that } \\
    & \forall z_0 \in D ~\&~ \forall n > n_0 ~\&~ \forall m \geq n ~\&~
              \forall z \in D\\
    & \left| v_{n,n}(z_0,z) - v_{m,m}(z_0,z) \right| < \epsilon.
  \end{align*}
\end{conj}

\subsection{Numerical Evidence}
\label{sec:evidence}

Figure with image of (0,0) in one step

Figures with areas corresponding to $V_n$ in \eqref{eq:volume} for two
possible images of (0,0) with $n=2$.

Figure of $V_2((0,0),z_b)$ vs $z_b$ indicated by color

Figure of $V_{5}((0,0),z_b)$ vs $z_b$ indicated by color

Figure with lines of $V_{5}((0,0),z_b)$ vs $z_b$ along the boundary
and along the diagonal that goes to the maximum on the boundary.

Repeat previous two figures for $V_{10}$

\section{Marginal Distributions}
\label{sec:marginal}

I want to describe the probability of sets of functions of a real
variable in terms of the stationary stochastic process corresponding
to \eqref{eq:phir}.  In this section I describe some marginal
probability distributions of that process.

\subsection{The Stationary Distribution $\stationary$}
\label{sec:convergence}

Does the stationary distribution $\stationary_d$ converge to a limit
as $d \rightarrow \infty$?
Recall \eqref{eq:stationary}
\begin{equation*}
  \stationary(z) = \eval \refunc(z) \lefunc(z),
\end{equation*}
and define $z_d \equiv (0,d)$.  The stationary distribution is
invariant under $\scale{}$ if and only if, for all $z$:
\begin{align}
  \frac{\stationary(z)}{\stationary(z_d)} &=
  \frac{\stationary\supn(z\supn)}{\stationary\supn(z\supn_d)} \nonumber\\
%
  \frac{\lefunc(z)\refunc(z)}{\lefunc(z_d)\refunc(z_d)} &=
  \frac{\lefunc\supn(z\supn)\refunc\supn(z\supn)}
       {\lefunc\supn(z_d\supn)\refunc\supn(z_d\supn)}
  \nonumber \\ 
  \label{eq:h_balance}
  \frac{\lefunc(z)\refunc(z)}
       {\lefunc\supn(z\supn)\refunc\supn(z\supn)} &= C.
\end{align}
Since $\lefunc((h,g)) = \refunc((-h,g))$, it is necessary that
\begin{equation}
  \label{eq:h_zero}
  \frac{\lefunc((0,g))}{\lefunc\supn((0,g\supn))} = \sqrt{C} ~\forall g
\end{equation}
hold for \eqref{eq:h_balance} to be true, and
\begin{equation}
  \label{eq:h_any}
  \frac{\lefunc(z)}{\lefunc\supn(z\supn)} = \sqrt{C} ~\forall z,
\end{equation}
which Conjecture \ref{conj} implies, is sufficient.

\subsection{Joint Distributions $\Skip{}$}
\label{sec:convergence_joint}

For a stationary stochastic process indexed by a continuous variable
$y$, the probability function for joint events\footnote{Strictly, I
  should write $\text{Prob}(z(y_1)\in A, z(y_2) \in B$ is a function
  of $y_2-y_1$.} is a function of their separation, ie,
\begin{equation}
  \label{eq:fdy}
  P(z(y_1)=z_a, z(y_2)=z_b) = f_{y_2-y_1}(z_a,z_b).
\end{equation}
If $y_2-y_1 = m\Delta_y$, one may approximate the function $f$ in
\eqref{eq:fdy} by $\Skip{m}(z_a,z_b)$.  After a rescaling,
$\scale{n}$, the appropriate approximation is
$\Skip{nm}\supn(z_a\supn,z_b\supn)$.  Combining \eqref{eq:skip},
\eqref{eq:h_any} and Conjecture~\ref{conj} yields
\begin{equation}
  \label{eq:flim}
  f_{y_2-y_1}(z_a,z_b) \equiv \lim_{n\rightarrow \infty}
  \lefunc\supn(z_a\supn) \pfop(y_2-y_1) \refunc\supn(z_b\supn).
\end{equation}

\section{Miscellaneous Properties}
\label{sec:properties}

\subsection{Boundary Points}
\label{sec:boundary}

Note that for $z_0 = (12,d-6)$,
\begin{equation*}
  F(z_0) = (0,d) \equiv z_d.
\end{equation*}
I assume that in a neighborhood of $z_d$ the eigenfunction is so
smooth that for any $z$ in a neighborhood of $z_0$ the approximation
\begin{equation*}
  \rho(z) = \frac{\rho(z_d)}{\lambda} a(z) \equiv \rho_0  a(z),
\end{equation*}
where $a(z)$ is the area of the image of $z$ is accurate to first
order.  In particular, for $z_a(\epsilon) \equiv (12-\epsilon, d-6)$
\begin{align*}
  z_{a1} &= (-\epsilon, d - 12 + (12-\epsilon)) \\
  &= (-\epsilon, d - \epsilon) \\
  z_{a2} &= (0,d) \\
  z_{a3} &= (\sqrt{24 \epsilon}, d - \epsilon) \\
  \rho(z_a(\epsilon)) &= \sqrt{6} \epsilon^{\frac{3}{2}}\rho_0.
\end{align*}
For $z_b(\epsilon) \equiv (12, d-6-\epsilon)$
\begin{align*}
  z_{b1} &= (0,d-\epsilon) \\
  z_{b2} &\approx (\epsilon ,d - \frac{\epsilon^2}{24}) \\
  z_{b3} &= (\sqrt{24 \epsilon}, d - \epsilon) \\
  \rho(z_b(\epsilon)) &= \sqrt{6} \epsilon^{\frac{3}{2}}\rho_0.
\end{align*}
And on the boundary with
\begin{align*}
  z_c & \equiv \left(12 - \epsilon , d - 6 + \epsilon -
  \frac{\epsilon^2}{24} \right) \\
  z_{c1} &= \left( -\epsilon, d - \frac{\epsilon^2}{24} \right) \\
  \rho(z_c(\epsilon)) & = \frac{\epsilon^3}{72} \rho_0.
\end{align*}
Note that since along the boundary at $z_0 = (12, d-6)$
\begin{equation*}
  \left. \frac{dg}{dh} \right|_{z_0} = -1,
\end{equation*}
the results for $z_a(\epsilon)$, $z_b(\epsilon)$ and $z_c(\epsilon)$
are consistent.

\subsection{Derivatives}
\label{sec:derivatives}

\begin{align*}
  \partialtwo{\rho(z_0)}{h_0} &= \frac{-1}{\lambda} \int_{h_1}^{h_3}
  \rho(h, g_1) dh \\
  \partialtwo{\rho(z_0)}{g_0} &= \frac{1}{\lambda} \left(
    \int_{h_1}^{h_3} \rho(h, g_1) dh
    - \int_{h_1}^{h_2} \rho(h, g_0 + 6 -h ) dh \right)
\end{align*}

\section{Speculation}
\label{sec:speculation}

In Figure~\ref{fig:eric} there is flow that is counter clockwise with
an overall drift to the right.  The approximation of the resulting
eigenfunction that appears in Figure~\ref{fig:Av2} varies by more than
a factor of $10^{40}$.  The eigenfunction is largest along boundary
where the flow accumulates which is on the right in
Figure~\ref{fig:eric}.  There seem to be the following three regimes:
\begin{description}
\item[Small $g$] Where images are truncated by the $g=-d$ boundary.
  Figure \ref{fig:Av2} shows that the truncation reduces the amplitude
  of eigenfunction along the $+h$ edge for small $g$.
\item[Intermediate $g$]
\item[Big $g$:] Near $g=d$ which is the only region where the flow
  maps points away from the right hand boundary.
\
\end{description}
I conjecture that for large $d$, the effect of the truncation at small
$g$ decays as the flow moves through the intermediate $g$ regime and:
\begin{itemize}
\item The eigenfunctions away from the small $g$ region converges as
  $d\rightarrow \infty$
\item The eigenvalue converges as $d\rightarrow \infty$ (Necessary for
  the eigenfunction convergence above.)
\end{itemize}
I don't know if the conjecture implies convergence to a measure on
functions as $\Delta_y \rightarrow 0$.

\section{Discrete Approximation}
\label{sec:approximate}

\subsection{2-d Illustrations}
\label{sec:ill}

I use $A$ to denote a discrete approximation of $\PFop$, and I define
$S$ by
\begin{equation*}
  S (h,g) \equiv (-h,g).
\end{equation*}
$A$ maps states $z_0 \equiv (g_0,h_0)$ to sets of states
$\left\{z_1 \right\}$, and while $A^T$ maps from $z_1$ to $z_0$ it is
not $A^{-1}$.  However,
\begin{align*}
  SS &= I \\
  SA^TS &= A,
\end{align*}
which ought to be worth something.

Figures \ref{fig:Av1}-\ref{fig:Av3} illustrate the action of $A$ and
$A^T$, and a numerically estimated eigenfunction appears in
Fig.~\ref{fig:eigenfunction}.  The figures depend on various values of
$\Delta_y$.  I should express them in coordinates
from~\eqref{eq:rescale}.
\begin{figure}
  \centering
    \resizebox{1.0\columnwidth}{!}{\includegraphics{Av1.pdf}}
    \caption{Image under $A^T$ of five points.}
  \label{fig:Av1}
\end{figure}

\begin{figure}
  \centering
    \resizebox{1.0\columnwidth}{!}{\includegraphics{Av2.pdf}}
    \caption{Image under $A$ of the same five points as used in
      Fig.~\ref{fig:Av1}.  For any point above the diagonal line in
      the upper left of the plot, the pre-image is truncated by the
      left edge of the allowed region.  That truncation explains the
      small values of the eigenfunction in that region.}
  \label{fig:Av2}
\end{figure}

\begin{figure}
  \centering
    \resizebox{1.0\columnwidth}{!}{\includegraphics{Av3.pdf}}
    \caption{Image of points under $A$ with a smaller value of
      $\Delta_y$.  Note that the the region of truncated pre-images is
      smaller than for Fig.~\ref{fig:Av2}.  However, since the point
      just to the right of that region has a pre-image that intersects
      the region, the eigenfunction will be suppressed there in a second
      order fashion.}
  \label{fig:Av3}
\end{figure}

\begin{figure}
  \centering
    %\resizebox{1.0\columnwidth}{!}{
      %\includegraphics{figs/2000_g_800_h_2_y_vec_log.png}}
    \caption{The $\text{log}_{10}$ of the amplitude of a numerically
      estimated eigenfunction.}
  \label{fig:eigenfunction}
\end{figure}
\end{document}

\onecolumn
\begin{landscape}
\section{Appendix}
\label{sec:appendix}

\newcommand{\rrrad}[3]{\sqrt{#1\left(#2-#3\right)}}
The following calculations demonstrate \eqref{eq:composition}:
\begin{equation*}
  \begin{matrix}
  \text{name} & \text{resolution }\Delta_y &\text{resolution }
  \frac{\Delta_y}{2} & \text{ name} & \text{showing}\\
  \hline
  z_0 & (h_0,g_0) & (2h_0, 4g_0) & z_0'=\scale{}(z_0)\\
  z_1=F(z_0) & (h_0-12,g_0+h_0-6) & (2h_0- 24, 4g_0 + 4h_0 -24) & z_1'
  = \scale{}(z_1)  \\
  && (2h_0-12, 4g_0+2h_0-6)& z_a'=F(z_0') \\
  && (2h_0 - 24, 4g_0 + 4h_0 -24) & F^2(z_0') & F^2 \circ \scale{} =
  \scale{} \circ F \\
  z_2 & (\rrrad{24}{d}{g_0} -12,g_0+\rrrad{24}{d}{g_0} -6)
  & (\rrrad{96}{d}{g_0} -24, 4g_0+\rrrad{192}{d}{g_0} -24) & z_2' =
  \scale{}(z_2) \\
  && (\rrrad{24}{4d}{4g_0} - 12, 4g_0 + \rrrad{24}{4d}{4g_0} - 6) &
  z_{a2}' \\
  && (\rrrad{96}{d}{g_0} - 24, 4g_0 + 2\rrrad{96}{d}{g_0} - 24 &
  F(z_{a2}') \\
  \scale{}^{-1}(F(z_{a2}')) & (\rrrad{24}{d}{g_0}
  -12,g_0+\rrrad{24}{d}{g_0} -6) &&&  \scale{}^{-1}(F(z_{a2}')) = \scale{}(z_2)
\end{matrix}
\end{equation*}
To do: Get sympy code in eric.py to calculate this and use results to make
an illustrative plot.  Perhaps include brief results from sympy here
rather than table above.
\end{landscape}

\twocolumn

\end{document}

first.py:    Code for eigenfunctions of integral operators for first
             order Markov process with 2-d states.

first_c.pyx: Cython version of LO from first.py

archive.py:  Initialize LO_step instance, calculate eigenfuction and
             store in archive directory

converge.py: Study convergence of numerical eigenfunction estimates

diff.py:     Plots the difference between two estimates of the Perron
             Frobenius function.  (Calculates the functions)

explore.py:  For looking at sensitivity of time and results to u, dy,
             n_g, n_h

map2d.py:    For making plots in 2d of regions to which points map and
             regions that map to points.

test.py:     Make sure that first.py and first_c.pyx make the same
             LO.  Uses archive

view.py:     Displays marginal densities and eigenfunctions using
             mayavi.  Uses archive
-----------------------------------------------------------------
Plan:

Symmetry is probably failing in s_bounds or ab.

Rethink grid.  Perhaps make (h=0, g=d) the reference point for the
grid

Get nosetest to call functions in test.py that ensure first and
first_c produce archive results and identical results to each other

Rewrite diff.py to use new LO_step class and write options for
comparing path integral calculations of relative volumes of paths from
a state to a pie slice.

nohup python3 conditional.py --d 16 --d_h .5 --d_g 1 --point 0 0
--iterations 1 2 4 8 16 32 48 --out test.pdf & 16 takes 51 hours, 32
dies

python3 conditional.py --d 16 --d_h 1 --d_g 2 --point 0 0 --iterations
1 2 4 8 16 # 16 takes 2.5 hours

nohup python3 conditional.py --d 16 --d_h 2 --d_g 4 --point 0 0
--iterations 1 2 4 8 16 32 & # 32 takes 7:19

%%%---------------
%%% Local Variables:
%%% eval: (TeX-PDF-mode)
%%% eval: (setq ispell-personal-dictionary "./localdict")
%%% End:
